program example_pstddm_conk_square
	use printVecMtx
	use PSTDDM
	use MySparseOperator
	! 	use mkl_pardiso
	implicit none
	integer, parameter :: sms = 4
	integer, parameter :: p = 2
	real(kind=8) :: eps = 1.d-10
	real(kind=8) :: O(2), width, height, hx, hy
	integer :: m, n, ms_pml(sms)
	complex(kind=8), dimension(:), allocatable :: u, ui, Fv

	real(kind=8), dimension(:), allocatable :: x, y
	integer :: ii, jj, err, im, m_pml, n_pml, kk
	
	real(kind=8) :: int_domain(4), sigma1, sigma2, k, pi
	integer :: m1, m2, n1, n2
	real(kind=8) :: errors(4)
	real(kind=8) :: dx_pml, dy_pml
	integer :: m_pml_seg, n_pml_seg, Mx, m_seg, Ny
	real(kind=8) :: dpar(4)
	integer :: ipar(6)

	integer :: dofs, nnz
	integer :: ios
	complex(kind=8), dimension(:), allocatable :: a
	integer, dimension(:), allocatable :: ia, ja
	real(kind=8) :: dy_seg, height_sub
	integer :: n_seg, n_sub
	
	pi = asin(1.d0)*2.d0

	O = (/ -0.24d+1, -0.24d+1 /)
	Mx = 10
	Ny = 10
	dx_pml = 0.4d0
	dy_pml = 0.4d0
	m_pml_seg = 1
	n_pml_seg = 1

	width = 2.d0*dx_pml + Mx*m_pml_seg*dx_pml
	height = 2.d0*dy_pml + Ny*n_pml_seg*dy_pml
	k = 0.314159265d+1 * 6	! wave number
	write(unit=*, fmt=*) 'k=', k
	
	sigma1 = 2.d+2 / k
	sigma2 = 2.d+2 / k
	int_domain = (/ O(1)+dx_pml, O(1)+width-dx_pml, O(2)+dy_pml, O(2)+height-dy_pml	 /)

	ms_pml = (/ 4, 8, 16, 32 /)

	err = 0

	do im = 1, sms, 1
		m_pml = ms_pml(im)
		n_pml = ms_pml(im)

		n = 2*n_pml + Ny*n_pml_seg*n_pml
		m = 2*m_pml + Mx*m_pml_seg*m_pml
		hx = width/m
		hy = height/n

		dpar(1:2) = O
		dpar(3) = dx_pml
		dpar(4) = dy_pml
		ipar(1) = m_pml
		ipar(2) = n_pml
		ipar(3) = m_pml_seg
		ipar(4) = n_pml_seg
		ipar(5) = Mx
		ipar(6) = Ny
		
		call pstddm_getdn(m,n,dofs)

! get the numerical solution without zeros boundary for the finite difference method 
		allocate(u(dofs), stat=err)
		if (err /= 0) print *, "u: Allocation request denied"


		allocate(Fv(dofs), stat=err)
		if (err /= 0) print *, "Fv: Allocation request denied"
		call pstddm_zrightvector(O,width,height,m,n, f_fun, Fv)
! get the solution of out PSTDDM
		call pstddm_conk_square(dpar,ipar,a_fun_full,c_fun_full,Fv,u)



		! get the interpolation of u_fun
! 		allocate(x(dofs), stat=err)
! 		if (err /= 0) print *, "x: Allocation request denied"
! 		allocate(y(dofs), stat=err)
! 		if (err /= 0) print *, "y: Allocation request denied"
! 		allocate(ui(dofs), stat=err)
! 		if (err /= 0) print *, "ui: Allocation request denied"
! 		forall(ii=1:m-1,jj=1:n-1)
! 			x(ii+(jj-1)*(m-1)) = O(1) + ii*hx
! 			y(ii+(jj-1)*(m-1)) = O(2) + jj*hy
! 		end forall
! 		call pstddm_zfdm(O,width,height,m,n, a_fun, c_fun, f_fun, ui)
! 		call u_funv(x,y,ui)
		
! get the errors of the numeircal solution for the FDM in the internal domain
		m1 = m_pml + 1
		m2 = m - m_pml - 1
		n1 = n_pml + 1
		n2 = n - n_pml - 1
		errors = (/	1.d0,	&	! get the L2 norm of u_fun
			&		1.d0,	&	! get the L2 error between u and  u_fun
			&		1.d0,	&	! get the H1 semi norm of u_fun
			&		1.d0 /)		! get the H1-semi errors bet
		call pstddm_errors(O,width,height,m,n,u_fun,uxy_fun,u,m1,m2,n1,n2,errors)
		call prtvec(errors,'f15.5')

		if (allocated(ui)) deallocate(ui, stat=err)
		if (err /= 0) print *, "ui: Deallocation request denied"
		if (allocated(y)) deallocate(y, stat=err)
		if (err /= 0) print *, "y: Deallocation request denied"
		if (allocated(x)) deallocate(x, stat=err)
		if (err /= 0) print *, "x: Deallocation request denied"
		if (allocated(u)) deallocate(u, stat=err)
		if (err /= 0) print *, "u: Deallocation request denied"
		if (allocated(Fv)) deallocate(Fv, stat=err)
		if (err /= 0) print *, "Fv: Deallocation request denied"
	end do


	
	if (allocated(Fv)) deallocate(Fv, stat=err)
	if (err /= 0) print *, "Fv: Deallocation request denied"
	if (allocated(y)) deallocate(y, stat=err)
	if (err /= 0) print *, "y: Deallocation request denied"
	if (allocated(x)) deallocate(x, stat=err)
	if (err /= 0) print *, "x: Deallocation request denied"
	if (allocated(u)) deallocate(u, stat=err)
	if (err /= 0) print *, "u: Deallocation request denied"
	if ( allocated(ui) ) deallocate(ui, stat=err) 
	if (err /= 0) print *, "ui: Deallocation request denied"
contains

	subroutine a_fun(x,y,a1,a2)
		real(kind=8), dimension(:,:), intent(in) :: x, y
		complex(kind=8), dimension(:,:), intent(out) :: a1, a2

		call a_fun_full(x,y,int_domain,a1,a2)
	end subroutine a_fun

	subroutine c_fun(x,y,c)
		real(kind=8), dimension(:,:), intent(in) :: x, y
		complex(kind=8), dimension(:,:), intent(out) :: c
		call c_fun_full(x,y,int_domain,c)
	end subroutine c_fun


	subroutine u_funv(x,y,u)
		real(kind=8), dimension(:), intent(in) :: x, y
		complex(kind=8), dimension(:), intent(out) :: u
		real(kind=8), dimension(:), allocatable :: r
		integer :: err,sx
		sx = size(x)
		allocate(r(sx), stat=err)
		if (err /= 0) print *, "r: Allocation request denied"
		r = dsqrt(x**2+y**2)
		where( r.le.eps )
			u = 0.d0
		end where
		where( r.le.1.d0 .and. r.gt.eps )
			u = r**3*(r*(r*(r+3.d0)-12.d0)+9.d0) * &
				&	dcmplx(bessel_j0(k*r), bessel_y0(k*r))
		end where
		where(r.gt.1.d0)
			u = dcmplx(bessel_j0(k*r), bessel_y0(k*r))
		end where
		if (allocated(r)) deallocate(r, stat=err)
		if (err /= 0) print *, "r: Deallocation request denied"
	end subroutine u_funv

	subroutine u_fun(x,y,u)
		real(kind=8), dimension(:,:), intent(in) :: x, y
		complex(kind=8), dimension(:,:), intent(out) :: u
		real(kind=8), dimension(:,:), allocatable :: r
		integer :: err,sx(2)
		sx = shape(x)
		allocate(r(sx(1),sx(2)), stat=err)
		if (err /= 0) print *, "r: Allocation request denied"
		r = dsqrt(x**2+y**2)
		where( r.le.eps )
			u = 0.d0
		end where
		where( r.le.1.d0 .and. r.gt.eps )
			u = r**3*(r*(r*(r+3.d0)-12.d0)+9.d0) * &
				&	dcmplx(bessel_j0(k*r), bessel_y0(k*r))
		end where
		where(r.gt.1.d0)
			u = dcmplx(bessel_j0(k*r), bessel_y0(k*r))
		end where
		if (allocated(r)) deallocate(r, stat=err)
		if (err /= 0) print *, "r: Deallocation request denied"
	end subroutine u_fun

	subroutine uxy_fun(x,y,ux,uy)
		real(kind=8), dimension(:,:), intent(in) :: x, y
		complex(kind=8), dimension(:,:), intent(out) :: ux, uy
		real(kind=8), dimension(:,:), allocatable :: r
		integer :: err,sx(2)

		if ( all(shape(x).eq.shape(ux)) ) then
			sx = shape(x)
			allocate(r(sx(1),sx(2)), stat=err)
			if (err /= 0) print *, "r: Allocation request denied"
			r = dsqrt(x**2+y**2)
			where ( r.lt.eps )
				ux = 0.d0
			end where
			where ( r.lt.1.d0 .and. r.ge.eps )
				ux = ( r**2 * ( 0.27d+2 + r*( -0.48d+2 + r*( 0.15d+2 + 0.6d+1*r)))*dcmplx(bessel_j0(k*r),bessel_y0(k*r))	&
					&	- k*r**3 * (0.9d+1 + r*( -0.12d+2 + r*( 0.3d+1 +r)))*dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) )	&
					&	* (-x/r)
			end where
			where ( r.ge.1.d0  )
				ux = k * dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) * x/r
			end where
		end if
		if ( all(shape(x).eq.shape(uy)) ) then
			if ( .not. allocated(r) ) then
				allocate(r(sx(1),sx(2)), stat=err)
				if (err /= 0) print *, "r: Allocation request denied"
				r = dsqrt(x**2+y**2)
			end if
			where ( r.lt.eps )
				uy = 0.d0
			end where
			where ( r.lt.1.d0 .and. r.ge.eps  )
				uy = ( r**2 * ( 0.27d+2 + r*( -0.48d+2 + r*( 0.15d+2 + 0.6d+1*r)))*dcmplx(bessel_j0(k*r),bessel_y0(k*r))	&
					&	- k*r**3 * (0.9d+1 + r*( -0.12d+2 + r*( 0.3d+1 +r)))*dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) )	&
					&	* (-y/r)
			end where
			where ( r.ge.1.d0  )
				uy = k * dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r)) * y/r
			end where
		end if
		ux = -ux
		uy = -uy
	end subroutine uxy_fun

	subroutine a_fun_full(x,y,int_domain,a1,a2)
		real(kind=8), dimension(:,:), intent(in) :: x, y
		real(kind=8), intent(in) :: int_domain(4)
		complex(kind=8), dimension(:,:), intent(out) :: a1, a2

		integer :: ii, jj, sx(2)

		sx = shape(x)

		do ii= 1, sx(1), 1
			do jj= 1, sx(2), 1
				if ( x(ii,jj).gt.int_domain(1) .and. x(ii,jj).lt.int_domain(2) .and. y(ii,jj).gt.int_domain(3) .and. y(ii,jj).lt.int_domain(4) ) then
					a1(ii,jj) = 1.d0
					a2(ii,jj) = 1.d0
				elseif ( x(ii,jj).le.int_domain(1) .and. y(ii,jj).gt.int_domain(3) .and. y(ii,jj).lt.int_domain(4) ) then
					a2(ii,jj) = dcmplx(1.d0,sigma1*(int_domain(1)-x(ii,jj))**p)
					a1(ii,jj) = 1.d0/a2(ii,jj)
				elseif ( x(ii,jj).ge.int_domain(2) .and. y(ii,jj).gt.int_domain(3) .and. y(ii,jj).lt.int_domain(4) ) then
					a2(ii,jj) = dcmplx(1.d0,sigma1*(x(ii,jj)-int_domain(2))**p)
					a1(ii,jj) = 1.d0/a2(ii,jj)
				elseif ( y(ii,jj).ge.int_domain(4) .and. x(ii,jj).gt.int_domain(1) .and. x(ii,jj).lt.int_domain(2) ) then
					a1(ii,jj) = dcmplx(1.d0,sigma2*(y(ii,jj)-int_domain(4))**p)
					a2(ii,jj) = 1.d0/a1(ii,jj)
				elseif ( y(ii,jj).le.int_domain(3) .and. x(ii,jj).gt.int_domain(1) .and. x(ii,jj).lt.int_domain(2) ) then
					a1(ii,jj) = dcmplx(1.d0,sigma2*(int_domain(3)-y(ii,jj))**p)
					a2(ii,jj) = 1.d0/a1(ii,jj)
				elseif ( x(ii,jj).le.int_domain(1) .and. y(ii,jj).le.int_domain(3) ) then
					a1(ii,jj) = dcmplx(1.d0,sigma2*(int_domain(3)-y(ii,jj))**p) / dcmplx(1.d0,sigma1*(int_domain(1)-x(ii,jj))**p)
					a2(ii,jj) = 1.d0/a1(ii,jj)
				elseif ( x(ii,jj).ge.int_domain(2) .and. y(ii,jj).le.int_domain(3) ) then
					a1(ii,jj) = dcmplx(1.d0,sigma2*(int_domain(3)-y(ii,jj))**p) / dcmplx(1.d0,sigma1*(x(ii,jj)-int_domain(2))**p)
					a2(ii,jj) = 1.d0/a1(ii,jj)
				elseif ( x(ii,jj).ge.int_domain(2) .and. y(ii,jj).ge.int_domain(4) ) then
					a1(ii,jj) = dcmplx(1.d0,sigma2*(y(ii,jj)-int_domain(4))**p) / dcmplx(1.d0,sigma1*(x(ii,jj)-int_domain(2))**p)
					a2(ii,jj) = 1.d0/a1(ii,jj)
				elseif ( x(ii,jj).le.int_domain(1) .and. y(ii,jj).ge.int_domain(4) ) then
					a1(ii,jj) = dcmplx(1.d0,sigma2*(y(ii,jj)-int_domain(4))**p) / dcmplx(1.d0,sigma1*(int_domain(1)-x(ii,jj))**p)
					a2(ii,jj) = 1.d0/a1(ii,jj)
				end if
			end do
		end do
	end subroutine a_fun_full

	subroutine c_fun_full(x,y,int_domain,c)
		real(kind=8), dimension(:,:), intent(in) :: x, y
		real(kind=8), intent(in) :: int_domain(4)
		complex(kind=8), dimension(:,:), intent(out) :: c

		integer :: ii, jj, sx(2)

		sx = shape(x)

		do ii= 1, sx(1), 1
			do jj= 1, sx(2), 1
				if ( x(ii,jj).gt.int_domain(1) .and. x(ii,jj).lt.int_domain(2) .and. y(ii,jj).gt.int_domain(3) .and. y(ii,jj).lt.int_domain(4) ) then
					c(ii,jj) = - k**2
				elseif ( x(ii,jj).le.int_domain(1) .and. y(ii,jj).gt.int_domain(3) .and. y(ii,jj).lt.int_domain(4) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma1*(int_domain(1)-x(ii,jj))**p)
				elseif ( x(ii,jj).ge.int_domain(2) .and. y(ii,jj).gt.int_domain(3) .and. y(ii,jj).lt.int_domain(4) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma1*(x(ii,jj)-int_domain(2))**p)
				elseif ( y(ii,jj).ge.int_domain(4) .and. x(ii,jj).gt.int_domain(1) .and. x(ii,jj).lt.int_domain(2) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma2*(y(ii,jj)-int_domain(4))**p)
				elseif ( y(ii,jj).le.int_domain(3) .and. x(ii,jj).gt.int_domain(1) .and. x(ii,jj).lt.int_domain(2) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma2*(int_domain(3)-y(ii,jj))**p)
				elseif ( x(ii,jj).le.int_domain(1) .and. y(ii,jj).le.int_domain(3) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma2*(int_domain(3)-y(ii,jj))**p) * dcmplx(1.d0,sigma1*(int_domain(1)-x(ii,jj))**p)
				elseif ( x(ii,jj).ge.int_domain(2) .and. y(ii,jj).le.int_domain(3) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma2*(int_domain(3)-y(ii,jj))**p) * dcmplx(1.d0,sigma1*(x(ii,jj)-int_domain(2))**p)
				elseif ( x(ii,jj).ge.int_domain(2) .and. y(ii,jj).ge.int_domain(4) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma2*(y(ii,jj)-int_domain(4))**p) * dcmplx(1.d0,sigma1*(x(ii,jj)-int_domain(2))**p)
				elseif ( x(ii,jj).le.int_domain(1) .and. y(ii,jj).ge.int_domain(4) ) then
					c(ii,jj) = (-k**2)*dcmplx(1.d0,sigma2*(y(ii,jj)-int_domain(4))**p) * dcmplx(1.d0,sigma1*(int_domain(1)-x(ii,jj))**p)
				end if
			end do
		end do
	end subroutine c_fun_full

	subroutine f_fun(x,y,f)
		real(kind=8), dimension(:), intent(in) :: x, y
		complex(kind=8), dimension(:), intent(out) :: f
		real(kind=8), dimension(:), allocatable :: r
		integer :: err
		allocate(r(size(x)), stat=err)
		if (err /= 0) print *, "r: Allocation request denied"
		r = dsqrt(x**2+y**2)
		where ( r.le.eps )
			f = 0.d0
		end where
		where(r.lt.1.d0 .and. r.gt.eps)
			f = -((((36.d0*r+75.d0)*r-192.d0)*r+81.d0)*r)*dcmplx(bessel_j0(k*r),bessel_y0(k*r)) &
				& + (((((12.d0*k)*r+(30.d0*k))*r-(96.d0*k))*r+(54.d0*k))*r**2)*dcmplx(bessel_jn(1,k*r),bessel_yn(1,k*r))
		end where
		where(r.ge.1.d0)
			f = 0.d0
		end where
		if (allocated(r)) deallocate(r, stat=err)
		if (err /= 0) print *, "r: Deallocation request denied"
	end subroutine f_fun


end program example_pstddm_conk_square










